## Loading required package: s20x

Handout 1: Introduction and null model

Handout 2: Simple linear model

Handout 3: Fitting curves using linear models

Handout 4: Categorical explanatory variable:

Handout 5: Multiplicative models — working with logged response variable 乘法模型

Handout 6: Categorical and numeric explanatory variables (a.k.a. ANCOVA model)

Handout 7: Categorical and numeric explanatory variables (a.k.a. ANCOVA model) 分类和数值变量

Handout 8: Models with several explanatory variables 多重回归模型

Handout 9: Power law models 幂律模型

Handout 10: 多级分类变量以及多级比较问题

Handout 11: 两个分类变量

Handout 12: 分析计数值数据

计数值数据

  • 计数值数据常出现在对分类数据的观察中
  • 如果对每个对象的度量都是分类变量(即:测量结果不是数值,而是归入那一类)
  • 分类的级级别就是categories(分类变量可以取到的不同值)
  • 计数值 count 就是每种分类组合出现的次数
  • rowdistr(.tbl) 查看行内百分比
  • count 数据不是正态分布的,线性回归有时候可以应用到count data(通常要log-transformed)

泊松分布

  1. non-negative 非负整数
  2. no upper limit 没有上限
  3. discrete variable 离散变量
  4. 重要特性 \(\sigma^2 = \mu\) : \(Var(Y)=E(Y)\)
  5. right skew when \(\mu\) is small, like normal when \(\mu\) is large 右偏态
    • 不看峰值,看尾部,尾部在哪边就是哪偏,右偏尾部为正直所以又称为正偏。
  6. cell counts 大于5的时候可以用卡方检测,将count当成正态分布而不是泊松分布
  • Useful functions
dpois(12, 9.61)
## [1] 0.08685078
rpois(20, 10)
##  [1] 10 17 12 10 12 10  6 11  9  6 13  9  6 12  8  7  9  5 15  7
barplot(dpois(1:12, 3), names=1:12, las=1)

barplot(dpois(0:250,100), names=0:250, las=1)

## 用卡方检测,检查2-way table的相关性 * 什么是卡方分布 若n个相互独立的随机变量ξ₁、ξ₂、……、ξn,均服从标准正态分布(也称独立同分布于标准正态分布),
则这n个服从标准正态分布的随机变量的平方和 \[Q = \sum_{i=1}^n \xi_i^2\] 构成一新的随机变量,其分布规律称为\(\chi^2\)分布 (原来某某分布指的就是随机变量值的分布规律)
(chi-square distribution),其中参数n称为自由度,正如正态分布中均值或方差不同就是另一个正态分布一样,自由度不同就是另一个\(\chi^2\)分布,记为\(Q \sim \chi^2(n)\) 或者 \(Q \sim \chi_n^2\).
卡方分布是由正态分布构造而成的一个新的分布,当自由度n很大时, \(\chi^2\)分布近似为正态分布。
对于任意正整数x, 自由度为 k的卡方分布是一个随机变量X的机率分布.

  • 期望和方差
    • 分布的均值\(\mu\)为自由度 n,记为 \(E(\chi^2) = n\)
    • 分布的方差\(\sigma^2\)为2倍的自由度(2n),记为 \(D(\chi^2) = 2n\)
    • 回忆:泊松分布的方差等于均值\(\sigma^2 = E(y) = \hat\mu\)
  • 性质
    • 分布在第一象限内(指的是概率密度函数),卡方值(随机变量每个值都是平方和)都是正值,呈正偏态(右偏态),随着参数 n 的增大,\(\chi^2\)分布趋近于正态分布;卡方分布密度曲线下的面积都是1.
    • 分布的均值与方差可以看出,随着自由度n的增大,\(\chi^2\)分布向正无穷方向延伸(因为均值n越来越大),分布曲线也越来越低阔(因为方差2n越来越大)。
    • 不同的自由度决定不同的卡方分布,自由度越小,分布越偏斜。
    • \(\chi_{v_1}^2\)\(\chi_{v2}^2\) 互相独立,则:\(\chi_{v_1}^2 + \chi_{v2}^2\) 服从\(\chi_{v_1+v_2}^2\)分布,自由度为 \(v_1+v_2\)

相关性检验逻辑:

  • 卡方检验\(H_0\)假设行列之间没有相关性;
  • 每一个Contingency Table里面的每个cell都有一个估值;
  • 一个subject处在该cell的几率:处在该行的几率×处在该列的几率估值 = \(\frac{n_{i+}}{n} \times \frac{n_{+j}}{n}\)
  • \(H_0\) 假设下这个cell的估值应该等于number of tatal counts(n) × 几率估值:\(\hat{n_{ij}} = n \times \frac{n_{i+}}{n} \times \frac{n_{+j}}{n}\)
  • 观察值预估值的差/噪音是否有足够的magnitude(量级)来证明\(H_0\)假设
  • 我们通过计算Z-statistic来测量H0
    • \(Z-statistic = \frac{n_{ij} - \hat{n_{ij}}}{sd(n_{ij})}\) 观测值-估值/观测值的标准差
    • 这里有一个原的假设,就是我们假设cells里的counts are Poisson distributed.
    • \(\Rightarrow sd({n_{ij}}) = \sigma = \sqrt{\mu}=\sqrt{E(Y)} = \sqrt{\hat{n_{ij}}}\)
    • \(Z_{ij} = \frac{n_{ij} - \hat{n_{ij}}}{\sqrt{\hat{n_{ij}}}}\) Z值is approximate符合正态分布 normal distribution
    • \[每个cell都有对应的Z_{ij}^2值,他们的和\sum_i\sum_j Z_{ij}\] 成为一个新的随机变量\(X\)
  • \(\hat{n_{ij}}\) 确定,\(Z-statistic\) 就确定, \(Z^2\) 确定, summing \(Z^2\) 得出 \(X\), 根据\(X\)和自由度求P-Value
    • 这个随机变量符合\(\chi_{(r-1)(c-1)}^2\) 的卡方分布 对于2-by-2 table 自由度为1
    • 由于pchisq(X,q) 可用来求随机变量\(\chi_q^2\)小于或等于分为点\(X\)的概率
    • 则:p-value = 1-pchisq(X,q) 对于2-by-2 table q=1
    • 也就是说 \[\sum_i \sum_j Z_{ij}^2\] 实际上是每个cell的估值-实际值相对偏差的累加效果的最佳评估,当其大到一定程度的时候我们推翻\(H_0\)假设,因为相关假设就是他们应该没有偏差。
    • 以上的每一步都可以被求出

有效条件:

  • 卡方检测用到了Z统计量,这就要求\(n_{ij}\) 正态分布:is at least approximately normally distributed.
  • \(\Leftarrow\) Expected counts \(\hat{n_{ij}}\) to be sufficient \(\geq 5\)

举例:

The \(\chi^2\) test of independence compares the oberved counts to thoes one would expect if there was no association between row and column factors/

AP.tb1 = matrix(c(17,83,27,19), nrow = 2, byrow = T, dimnames = list(c("attand", "not.attend"),c("fail","pass")))
AP.tb1
##            fail pass
## attand       17   83
## not.attend   27   19
rowdistr(AP.tb1)
## Row Proportions
##            fail pass Totals   n
## attand     0.17 0.83      1 100
## not.attend 0.59 0.41      1  46

# show to show row proportions
chisq.test(AP.tb1)
## 
##  Pearson's Chi-squared test with Yates' continuity correction
## 
## data:  AP.tb1
## X-squared = 24.073, df = 1, p-value = 9.274e-07
chisq.test(AP.tb1, correct = F)
## 
##  Pearson's Chi-squared test
## 
## data:  AP.tb1
## X-squared = 26.016, df = 1, p-value = 3.386e-07

默认chisq.test 是修正的,这个修正是为了通过修正\(Z_{ij}\)的非连续性,来修正Counts的非连续整数的性质
使用卡方检验的条件是:count不能太小,太小就不正态分布了,\(n_{ij}>5\) 如果小于则代码会有warning 1. P-Value表明了Null hypothesis was rejected. 他们之间有相关性

AP.chisq=chisq.test(AP.tb1, correct = F)
AP.chisq$expected
##                fail     pass
## attand     30.13699 69.86301
## not.attend 13.86301 32.13699
X = (17-30.14)^2/30.14 + (83-69.86)^2/69.86 + (27-13.86)^2/13.86 + (19-32.14)^2/32.14
X
## [1] 26.02961
chisq.test(AP.tb1, correct = F)
## 
##  Pearson's Chi-squared test
## 
## data:  AP.tb1
## X-squared = 26.016, df = 1, p-value = 3.386e-07
  1. 从估值也能看出他们的\(H_0\)假设不成立

The next step is to quantify the strength of the association.

rowdistr(AP.tb1, comp = "between", plot = F)
## Row Proportions
##            fail pass Totals   n
## attand     0.17 0.83      1 100
## not.attend 0.59 0.41      1  46
## 
## 95% CIs for diffs between proportions with fac2 = fail 
## (rowname-colname) 
##        not.attend      
## attand (-0.577,-0.257)*
## 
## 95% CIs for diffs between proportions with fac2 = pass 
## (rowname-colname) 
##        not.attend    
## attand (0.257,0.577)*
rowdistr(AP.tb1, comp = "between", suppressText = T)

通过比较同一列在不同行内的概率的不同,可以知道行列有相关性(interactionPlot 不行) 但并不能比较相关性的强弱,即概率差相同但在尺度上不同,不能认为他们的相关度就相同

卡方检验的问题:

  1. count < 5 会报错,因为cell里的count
  2. 需要chisq.test(*.tb, correct=F), TRUE和F要对比着看,如果差异太大也有问题 默认是True,因为需要修正Zij的非连续性,以便修正count不是连续数的问题
  3. chisq.test()$expected 可以查看估值
  4. rowdistr(*.tb, comp=“between”, plot=F) 可以查看同一列内不同的概率,以及95% CI 概率差估值,根据概率不同,可以看到他们的有关联度。
  5. 根据概率差的大小可以看出他们的关联度强弱
  6. 同一列内不同行的概率差的置信区间(xxx, xxx)这样的区间在2-by-2 contingency table里有两个,在column为two levels factor 时他们绝对值相同互为正负
  7. 这个置信区间可以手动计算
    • \(se(\hat p_1 - \hat p_2) = \sqrt{\frac{\hat p_1(1-\hat p_1)}{1} + \frac{\hat p_2(1- \hat p_2)}{1}}\)
    • \(\hat p_1 - \hat p_2 \pm 1.96 \times se(\hat p_1 - \hat p_2)\)
    • 这个概率差是总结时需要用到的。
  8. 这个不同row的概率差,在不同的column有所不同,提供了evidence of association. 但其很难在尺度上衡量association的strength。
  9. interactionPlots()不同,他的值是统计次数,还要看斜率才能粗略知道是否相关。

Handout 13: 比值比

AP.tb1
##            fail pass
## attand       17   83
## not.attend   27   19
  1. 手动求CI for OR的公式 知道contingency table 可以交叉求出\(\hat{OR}\),
    注意排列顺序不同\(\hat{OR}\)可能不同,但也可能相同。
  2. \(\hat{OR} = \frac{\hat{Odds_1}}{\hat{Odds_2}} = \frac{n_{11}n_{22}}{n_{12}n_{21}}\)
    • 这个\(\hat{OR}\) 算出来就可以知道同在列1组内不同行事件概率的比值了,而且可以推广到不同table排列方式的比值
    • \(se\left(\log(\hat{OR})\right) = \sqrt{\frac{1}{n_{11}} + \frac{1}{n_{12}} + \frac{1}{n_{21}} + \frac{1}{n_{22}}}\) 根据STATS703
    • 当样本比较大的时候\(\log(\hat{OR})\) 接近正态分布 is approximately normal,
    • so CI for \(\log(\hat{OR})\):
      • \(\Rightarrow \log(\hat{OR}) \pm 1.96 \times se\left(\log(\hat{OR})\right) = (a,b)\)
      • \(\Rightarrow\) 95% CI for \(\hat{OR}\):
      • \(\exp(a), \exp(b)\) 用以解释数据的行列比较
  3. Test the association using the log odds: \(H_0\) means no row-column association \(\Rightarrow \hat{OR}=1 \Rightarrow log(\hat{OR})=0\) To test \(H_0\) we use Z statisic \(Z = \frac{\log(\hat{OR}) - 0}{\rm se(\log(\hat{OR}))}\) \(\Rightarrow P-value 进而判断 \hat{OR}=1 是否极端\) 其实这个P-Value就是测试\(H_0:log(\hat{OR}=0)\) 跟第15章的 summary(glm.fit) 出来的第四二个 \(P-value 测试 coefficient[4] \beta_3=0\) 是一样的,p值也一样 z-test 适用于大型数据,多余30个个体以上,如果知道标准差,那最好用z!

Handout 14: 费舍尔检验 FET

Handout 15: 推广线性统计模型到关联列表

  1. fisher.test(tbl)$p.val \(H_0\) 是无相关 或者 `chisq.test(tbl)$expected 倍数差距相等 则无相关
  2. 模型: \(n_{11} = n \times r_1 \times c_1\)
    \(log(n_{11}) = log\,n + log\, r_1 + log\, c_1 = \beta_0\)
    \(log(n_{12}) = \beta_0 + \beta_1 \Rightarrow \beta_1 = log c_2 - log c_1\)
    \(log(n_{21}) = \beta_0 + \beta_2 \Rightarrow \beta_2 = log r_2 - log r_1\)
    \(log(n_{22}) = \beta_0 + \beta_1 + \beta_2\) no interaction
    \(log(n_{22}) = \beta_0 + \beta_1 + \beta_2 + \beta_3\) interaction
    \(\beta_3 = log(\hat{OR})\)
  3. fit = glm(y~fac1 * fac2, family = poisson, data = df)
  4. summary(fit) residual deviance means residual sum-of-squares 4. summary 出来就coefficient \(\beta_0, \beta_1, \beta_2\) 将他们加起来再exp 就可以得到4个格子的期望值了。似乎没什么意思,估出来跟观测值相同,因为模型中我们加入了相关性参数。
  5. 如果\(\beta_3 \neq 0\)
    • 其中\(\beta_3\) 的解释比较多 取exp 可以得到 \(\hat{OR}\)
    • 用起来解释不同组合的odds之比
  6. 如果\(\beta_3 = 0\)
    • 没有interaction的情况下,我们可以说,独立的针对行内和列内进行比较
    • \(log(n_{12}) = \beta_0 + \beta_1 \Rightarrow \beta_1 = log c_2 - log c_1 \Rightarrow exp^{\beta_1} = \frac{c_2}{c_1}\)
    • \(log(n_{21}) = \beta_0 + \beta_2 \Rightarrow \beta_2 = log r_2 - log r_1 \Rightarrow exp^{\beta_2} = \frac{r_2}{r_1}\)
    • c2列出现的数量是c1列出现数量的多少倍
    • r2行出现的数量是r1行出现数量的多少倍

Handout 16: 推广线性统计模型到计数数据

举例1:

  • 地震的例子其实是:一个 ANCOVA with interaction
  • 对于glm(fac1*x2, family = poisson)
  1. 由于\(log(E[Y]) 而且 \beta_2 是相对于 reference 的斜率\)
  2. 要求另一组斜率就要调换fac1 里面两个levels的顺序,或者\(\beta_2 + \beta_3\)得出另一个斜率

举例2:

  • 大学录取性别歧视检测其实是一个 two-way ANOVA
  • 最开始的时候只考虑性别和是否录取:two-way ANOVA
  • 只不过原始数据没有给出关联列表,而是要用xtabs(y~gender+outcome) 得出关联列表
  • 后来又考虑到学院就开始要用xtabs(Y~x1 + x2 + x3, data=df) 了
    • 这里有个著名的辛普森悖论,即,不同组里出现的趋势,在合并这些组以后改变了。
    • 模型变成了3个解释变量,变成\(fac_1 \times x_3 + fac_2 \times x_3 + fac_1 \times fac_2\)

举例3:

  • 鲷鱼在两个地区及在其保护区与否,的不同数量count
  • interactionPlots(y~fac_1 \times fac_2) 斜率相等不能证明没有相关性了。要在倍数尺度(multiplicative scale)上相等才行
  • \(glm(y~ fac_1 * fac_2)\)
  • summary(glm) \(residual deviance \approx df \Rightarrow no need quasi+\)
  • \(P-value < 0.05 \Rightarrow no association\)
  • Explain:
    • \(\beta_1: c2/c1\)
    • \(\beta_2: r2/r1\)

Handout 17.5 逻辑回归

当响应变量是二元量的时候 -(逻辑回归logistic regression)

  • Y是伯努利随机变量,它的取值为1/0,该伯努利实验有个参数:成功率(Y=1)的几率为p,
  • \(E(Y)=p\) 伯努利随机变量的均值 = p
  • 为了fit model:如何将区间在[0,1]的二元数据扩展到为\([-\infty, +\infty]\)
  • \(Odds = \frac{p}{1-p}\) \(\subset [0, +\infty]\)
  • log-odds is \(\log(\frac{p}{1-p}) = x\) \(\subset [-\infty, +\infty]\) 为logit:分对数函数
  • \(\Rightarrow p = \frac{exp(x)}{1+exp(x)}\) 为logistic

  • 技术上来讲,它也就是一个二项式的广义线性模型 二项式广义线性模型也可以用在概率上,只要你知道尝试的总数。 成功次数加起来除以总数就得到概率
  • 逻辑函数
    • 具体步骤:
    • 用xtabs()帮你统计成功数,并用:100*成功数/次数=得出成功率,但这不重要(每个性别在每个距离上都是投了10次,所以trial 次数为10)
    • \(logit(E[Y_i]) = \beta_0 + \beta_1 x_1 + \beta_2 x_2 + \beta_3(x_1 \times x_2)\)
    • fit = glm(Y~X_1 \times X_2, family = binomial, data = df)
    • \(Y_i\)是值为0,1的随机变量 通过取它的odds的对数而不是直接取它的对数,变为正态分布
    • 书本的男女不同距离投篮例子其实还是一个 two-way ANOVA
    • anova 这里给出的信息非常多:
    • p-value >0.05 的系数可以被简化掉了。基本上这样就够了
    • 1-pchisq(Resid.Dev (Residual Deviance) Resid.Df (Residual Df)) 小于.05 则用 quasi+
    • 根据这两个数可以求p-value = 1-pchisq(RD, DF) 当做泊松分布来计算行列相关性,>0.05则说明没有相关性得到的结果应该是很接近的
    • 也可以根据DF 求卡方分布的标准差\(s = \sqrt{2*DF}\)
  • refit的model
  • summary(fit, family=binary, data =)
  • 将系数还原\(exp(\beta)\) 得到的值是对Odds的影响
  • 分析个系数的意义 So a 1 metre increase in distance would result in a 92.8% reduction in odds. Base on the p-value from the summary, simplify the fitted model. exponentitate the coefficient to get the multiplicative effect on the odds

当相应变量是成功数的时候

  • 用xtabs(y~x, df = ) 统计成功数
  • 只是因变量不是[0,1]指,而是1的个数:the count of ‘successes’
  • 此处将因变量变化为:\(fit=glm(cbind(success, fail)~x, family = binormial, data = df)\)
  • 其他相同
  • 还是一个逻辑回归 logistic regression

Handout 18: 当相应变量是百分比的时候

Handout 18.5: 垃圾邮件

Handout 时间序列

介绍时间序列

时间序列分析:

  • 一个时间序列由在连续递增的时间上对数据的收集、重新排序和观察组成。
  • 时间间隔相等:时、日、月、季度、周年等
  • 时间序列的数据是有序的,不像回归,这里数据的顺序很关键
  • 注意,时间是有方向的,如果过去能够影响将来,那么观察值就不是相互独立的。

研究方法

  • head(data.df, 5)/tail 查看变量、值的类型
  • pairs20x(data.df[,2:5]) 为数据帧绘制直方图(对角线)和散点图(右上角)
  • 没有把时间放进去研究的时候
  • plot(data.df, which = 1) 查看残差分布,constant and with no trend
  • normcheck(lm.fit)
  • cooks20x(lm.fit) no influential outlier
  • 由于曲线弯曲,所以log是试试
  • refit.fit = lm(log(y)~x, data.df)
  • summary(lm.fit)
  • 查看残差随自变量的分布趋势
  • 把时间放进去一起研究的时候残差图看起来更像随机分布在0附近了。可见年顺序在时间序列分析中的作用。

时间序列的数据结构

  • 随着时间记录下来的响应变量值
  • 以及一系列在同样时间点记录下来的解释变量(表示为:\(X_{t1},X_{t2},X_{t3},...,X_{tk},\))
    • 我们将用它们解释相应变量的行为,或者预测相应变量未来的值。
  • 在金融和经济领域,我们通常关心股票市场、消费价格等随时间的变化。
  • 这类变量的观察值常常以指数形式出现:
    • 消费价格指数 CPI
    • 交易权重指数 TWI
    • 股指 NZSX50
  • 附带解释变量的观察值
    \(时间段 \quad 观察值 \quad 解释变量\)
    \(\quad 1 \qquad Y_1 \qquad X_{11} \quad X_{12} \quad \cdots \quad X_{1k}\)
    \(\quad 2 \qquad Y_2 \qquad X_{21} \quad X_{22} \quad \cdots \quad X_{2k}\)
    \(\quad \vdots \qquad \vdots \qquad \quad \vdots \qquad \vdots \quad \; \vdots \qquad \vdots\)
    \(\quad t \qquad Y_t \; \qquad X_{t1} \; \quad X_{t2} \quad \cdots \quad X_{tk}\)
    \(\quad \vdots \qquad \vdots \qquad \quad \vdots \qquad \vdots \quad \; \vdots \qquad \vdots\)
    \(\quad T \qquad Y_T \qquad X_{T1} \quad X_{T2} \quad \cdots \quad X_{Tk}\)
  • 这里多了一列时间,它的顺序必须是升序排列

时间序列的成分

时间序列模型

在为时间序列数据建模的时候,分解成四个成分可以帮助理解和量化数据随时间的变化 例如:我们的模型可以如下:
\[Y_t = T_t + C_t + S_T + R_t\]

  • 我们可以将这几个成分重新合并成一个新模型来做预测,但这不是理想的方法。

四个成分

  1. 趋势成分:
    最慢的成分,我们的问题是:这个趋势是上升的还是下降的?
    • 最慢、长期、平滑
    • 趋势通常来自总体的增长、变化、技术性的改变等等。
    • 趋势可以随着时间或者在一段时间内,上升或下降。
  2. 周期成分:
    不规则波浪形行的模式可以出现在不同的时间间隔,并且每次出现都不一样。
    • 第二慢,周期性的,一系列波浪形的
    • 波动的波长和振幅既不固定也无法预测,所以是不规则的模式
    • 很多时间序列的周期性行为主要来自商业周期 注:如果趋势成分允许非线性,则很难区分趋势成分和周期(波动)成分,这样他们就可以合并成一个趋势成分。
  3. 季节成分:
    有规律的重复的模式,有固定长度的时间间隔,例如周年或者月。
    • 变化比周期更迅速,为时间序列中,相等间隔规律性的重复模式。
    • 季节性基于日历或者小时
    • 可以字面上理解为三个月(一季度),也可以理解为月,周,日或年。
  4. 随机成分(余数)
    无法预测的随机变化
    • 总是存在的,变化最快的成分。
    • 它由那些无法预测的关于趋势,周期,季节成分的变化组成。
    • 但他可以当作其他模型的观测值一样拟合,例如标准正态分布\(N(0, \sigma^2)\), 但它可以不满足独立性原则。
    • 它可以理解为当其他成分都去掉以后留下的成分
    • 它是时间序列中无法被其他东西解释的成分
    • 在季节变化的数据中,随机成分可以从季节成分无法准确重复出现看出。

STL(Seasonal Trend Loess)函数分解时间序列

  1. 将响应变量转换成时间序列类型对象 ts()
  2. stl()函数执行季节趋势回归分解
data("airpass.df")
passengers.ts = ts(airpass.df$passengers, start = 1949, frequency = 12)
plot(passengers.ts, main = "", ylab = "")

passengers.stl = stl(passengers.ts, s.window = "periodic")
plot(passengers.stl)
abline(v=1950:1960, lty=2)

如何读图

  • 解释变量和它的值取自连续的时间间隔
    • 响应变量相对于时间的图有助于发现是否存在时间趋势。
plot(passengers.ts, main = "", ylab = "")

remainder = passengers.stl$time.series[,3]
    • 有对时间进行回归分析得出的残差图可以显示出其他解释变量所没有捕捉到的时间的作用。
    • 将观察值的点连接起来有利于发现时间顺序的模式。plot(y~x, data=data.df) lines(y~x, data=data.df)

例1:

data("mening.df")
mening.ts = ts(mening.df$mening, frequency=12, start=c(1990,1))
plot(mening.ts, main="")
abline(v=1991:2001, lty = 2)

# airpass.df = within(airpass.df, {month = factor(rep(month.abb, 12), levels = month.abb); time = 1:144})
# lines(passengers~month, data = airpass.df)

解读:

  • 可以看出有年度季节性模式
  • 每隔12个月放一个垂直线,将数据按成年分段:abline(v=1991:2001, lty=2)
  • 总体来看这个图是在重复的。
  • 所以有季节性。

例2

data(beer.df)
beer.ts = ts(beer.df, frequency = 12, start = c(1971,7))
plot(beer.ts, main = "US Beer Production, July 1970 to June 1978")
abline(v = 1971:1978, lty = 2)

解读:

  • 上图总体来看也是重复波动的,而且它的常规峰值出现在一年中相同的时间

例3

passengers.ts = ts(airpass.df$passengers, start = 1949, frequency = 12)
plot(passengers.ts, main = "", ylab="")
abline(v=1950:1960, lty=2)

解读:

  • 我们可以看到季节性的涨落幅度更大了;
  • 如果这个幅度随时间增加或者减少,我们应该用乘法分解代替加法分解。 \(Y_t = T_t \times C_t \times S_t \times R_t\)
  • 由于响应变量是一个计数值,所以用乘法模型是很自然的
  • 原则上,我可以通过扩展处理泊松数据的方法,来将其用在计数值时间序列上。
  • 这里我们用乘法模型,通过\(log(Y_t)\)来实现,前提是\(Y_t\)都比较大,如本例中。
  • 对数据取自然对数可以是季节变化变得恒定,但时间成分的作用将变成倍数(增加多少时间->数据增加多少倍)
log.passengers.ts = log(passengers.ts)
plot(log.passengers.ts, main = "", ylab="")
abline(v=1950:1960, lty=2)

分解图:

  • 分解图将帮助我们去掉季节和随机成分,去了解潜在的趋势或波动:
passengers.stl = stl(log.passengers.ts, s.window = "periodic")
plot(passengers.stl)
title("Decomposition of log(Passengers)")

* 趋势存在少量的非线性,最多表明数据中可能有较弱的波动。

分解时间序列:

要解决的问题:

  • 相对12月,1月的销量往往有所下降,我们想知道这个下降是否仅仅是正常的季节变化引起的?
  • 同样的,由于学生毕业加入劳动力人群,在每个学年末失业率都有所提升。年末提升的失业率表明失业率真的有所提升还是仅仅是正常的季节变化?

解决的方法:

  1. 研究去季节化后的值(时间序列去掉季节成分)
    1. 我们估计去季节后的值存在趋势或波动(或者存在解释变量的作用)。研究潜在的趋势或者波动
    2. 如果没有证据证明趋势的存在,我们就可以推断我们看到的只是时间序列上正常的季节变化
    3. 否则我们就说存在根本性的变化,而不是季节
  2. 用stl分解
    我们从数据中减掉季节成分再在其上做回归分析;注:这里要小心,因为不满足独立假设。
    1. 获得季节成分
    2. decomp.stl = stl(data.ts, s.window = “periodic”)
    3. sa.data = data.ts - decomp.stl$time.series[,2]
# 代码段2 ----
passengers.ts = ts(log(airpass.df$passengers), start = 1949, frequency = 12)
plot(passengers.ts, main = "Data with Seasonal Component")

passengers.stl = stl(passengers.ts, s.window = "periodic")
sa.passengers = passengers.ts - passengers.stl$time.series[, "seasonal"]
plot(sa.passengers, main = "Data with seasonal component removed")

用一个解释分类变量来去掉季节成分

  • 将一个分类变量用在时间序列的回归模型中,分类的级别对应每个季节。
  • 用这种方法不需要调整数据去掉季节成分,季节成分是模型设定的一部分
  • 这种回归方法适合用在预测,因为它考虑到了整个时间序列。
  • 另一种预测方法是Holt-Winters 指数平滑。

预测

指数平滑(exponential smoothing)

概念

  1. 指数平滑是能用过去的所有信息来预测将来的方法
  2. 每个平滑值是序列中之前所有值的加权平均数(weighted average)
  3. 指数序列标记为:\(S_t\) (区别于之前的原始序列\(Y_t\))
  4. \(S_t\) 的平滑度有平滑参数\(\alpha\)控制

基本形式:

\[S_t = \alpha Y_t + (1 - \alpha) S_{t-1}\]

  • 平滑常量\(\alpha\)越小,结果序列越平滑
  • 之所以叫指数平滑是因为距离现在越远的观察值其权重随时间指数级的变小(exponentially smaller)
  • 表示如下: \[S_t = \alpha Y_t + \alpha (1 - \alpha) Y_{t-1} + \alpha (1 - \alpha)^2 Y_{t-2} + \cdots + (1-\alpha)^{t-1} Y_1\]
  • 从等式的左往右,随着t变得非常大,表达式的系数变得更小(以指数速度 at an exponential rate)

性质

  • 平滑后的序列有赖于之前的所有值,其中离现在最近的值的权重最大
  • 指数平滑需要大量的观测值
  • 以上指数平滑的基本形式不适合存在趋势或者季节性的数据
  • 指数平滑在R中用 HoltWinters 函数

指数平滑的基本形式:(废除其他成分,免得成为通用形式)

  • HoltWinters(data.df, alpha = x, beta = FALSE, gamma = FALSE) \(\alpha\) 由x来制定,beta 和 gamma 用在有趋势和季节性的时候
data("larain.df")
LArain.ts = ts(rev(larain.df$LA.Rain) * 25.4, frequency = 1, start = 1877)
plot(LArain.ts, xlab = "Season", ylab = "Total rainfall (mm)")

LArain.hw = HoltWinters(LArain.ts, alpha = 0.5, beta = FALSE, gamma = FALSE)
plot(LArain.hw, lwd = 2, col.predicted = "blue", main = "alpha = 0.5")

LArain.hw = HoltWinters(LArain.ts, alpha = 0.95, beta = FALSE, gamma = FALSE)
plot(LArain.hw, lwd = 2, col.predicted = "blue", main = "alpha = 0.95")

LArain.hw = HoltWinters(LArain.ts, alpha = 0.75, beta = FALSE, gamma = FALSE)
plot(LArain.hw, lwd = 2, col.predicted = "blue", main = "alpha = 0.75")

LArain.hw = HoltWinters(LArain.ts, alpha = 0.25, beta = FALSE, gamma = FALSE)
plot(LArain.hw, lwd = 2, col.predicted = "blue", main = "alpha = 0.25")

LArain.hw = HoltWinters(LArain.ts, alpha = 0.05, beta = FALSE, gamma = FALSE)
plot(LArain.hw, lwd = 2, col.predicted = "blue", main = "alpha = 0.05")

  • 实际上我们可以让HoltWinters来指定\(\alpha\),它通过最小化均方预测误差来确定。

预测举例:

predict 中调用 n.ahead 参数

LArain.hw1 = HoltWinters(LArain.ts, beta = FALSE, gamma = FALSE)
LArain.pred = predict(LArain.hw1, n.ahead = 5)
plot(LArain.hw1, predicted.values = LArain.pred, lwd = 2, col.predicted = "blue")

LArain.pred = predict(LArain.hw1, n.ahead = 15, prediction.interval = TRUE)
plot(LArain.hw1, predicted.values = LArain.pred, lwd = 2, col.predicted = "blue")

hist(resid(LArain.hw1))

  • 这个模型似乎不是很好,对数据取对数来去掉向右歪斜,同时避免第
log.LArain.ts = log(LArain.ts)
LArain.hw2 = HoltWinters(log.LArain.ts, beta = FALSE, gamma = FALSE)
plot(LArain.hw2)

log.LArain.pred = predict(LArain.hw2, n.ahead = 15, prediction.interval = TRUE)
exp(log.LArain.pred)
## Time Series:
## Start = 1943 
## End = 1957 
## Frequency = 1 
##           fit      upr      lwr
## 1943 343.7816 929.1857 127.1928
## 1944 343.7816 931.5911 126.8644
## 1945 343.7816 933.9966 126.5377
## 1946 343.7816 936.4020 126.2126
## 1947 343.7816 938.8074 125.8893
## 1948 343.7816 941.2129 125.5675
## 1949 343.7816 943.6184 125.2474
## 1950 343.7816 946.0239 124.9289
## 1951 343.7816 948.4296 124.6121
## 1952 343.7816 950.8353 124.2968
## 1953 343.7816 953.2411 123.9831
## 1954 343.7816 955.6470 123.6709
## 1955 343.7816 958.0531 123.3604
## 1956 343.7816 960.4593 123.0513
## 1957 343.7816 962.8657 122.7438

用HoltWinters预测

  • 可以用HoltWinters预测包含趋势和季节的时间序列
  • 在拟合过程中有三个不同的指数平滑,更新公式:
    • 水平 作用 \(a_t\)
    • 斜率 作用 \(b_t\)
    • 季节 作用 \(s_t\)
  • 每一个更新的方程都有它自己平滑常数:\(\alpha, \beta , \gamma\)
  • 预测模型为: \[F_{T+h} = a_T + hb_T + s_{T,h}\]
    • 预测是通过估计水平,斜率和季节在时间上的作用来实现的。T是时间序列中的最后一个时间点。
    • \(s_{T,h}\) 表示持续重复的季节作用,尤其是在时间T上估算出来的,并用来计算h个时间单位以后的季节作用。
log.pass.ts = ts(log(airpass.df$passengers), frequency = 12, start = 1949)
log.pass.hw = HoltWinters(log.pass.ts)
log.pass.pred = predict(log.pass.hw, n.ahead = 30)
plot(log.pass.hw, log.pass.pred, col.predicted = "blue", lwd = 2)

让R来决定最优的 \(\alpha, \beta , \gamma\)

log.pass.pred = predict(log.pass.hw, n.ahead = 30, prediction.interval = TRUE)
plot(log.pass.hw, log.pass.pred, col.predicted = "blue", col.intervals = "red", lwd = 2)

针对时间序列的线性回归模型

时间序列线性回归

  • 我们前面一直假设在标准模型中,误差是相互独立的,但这个假设在时间序列不成立。
  • 通过明确建立依赖模型,可以扩展线性模型来放宽这个假设要求
  • 在此简述一阶(first-order)自相关作用

自相关作用

  • 拟合线性模型的到时间序列的主要问题就是随机误成分的值常常非独立
  • 这个问题成为自相关或者连续相关
  • 自相关很难人言识别,尤其是有很多解释变量的时候。所以要集中注意力在解释那些拟合解释变量到响应变量以后所剩下的产值上。

识别自相关

最简单的就是检查残差的自相关函数图(ACF) 首先要明确用回归模型分析时间序列的策略:

  1. 以可用的解释变量给均值建模,如果没有,到第二步
  2. 通过增加一个季节分类变量来给季节性(如有)建模(这个方法的缺陷是它没有施加任何模式到季节成分上,所以它没有利用1,2月比1月和7月更近这样的事实)
  3. 通过增加一个(或多个)时间项到模型中来给所有时间依赖关系建模
  4. 检查残差的ACF图。如果有自相关证据,延迟相应(lagged-response)模型可以解决

举例:

plot(log.passengers.ts)

  1. 没有解释变量,跳过第一步
  2. 有很强的季节效果,建立一个以季节分类建立一个线性模型
# 第二步:添加季节项 ----
airpass.df = within(airpass.df, { month = factor(rep(month.abb, 12), levels = month.abb); time = 1:144})
pass.fit = lm(log(airpass.df$passengers)~month, data = airpass.df)
plot(pass.fit$residuals)
# plot(pass.fit$residuals, type = "l")
lines(airpass.df$time, residuals(pass.fit)) # 用线将点连接起来
lines(lowess(airpass.df$time, residuals(pass.fit))) # 用离散点平滑后画出的曲线,用来了解趋势再好不过了。

# lowess returns a list containing components x and y which give the coordinates of the smooth. The smooth can be added to a plot of the original points with the function lines: see the examples. 在这里显然不用了,因为散点图已经可以看出明显的趋势了。
  1. 可以看到很强的时间作用(线性),需要再对时间建模,
# 第三部:添加时间项 ----
pass.fit1 = lm( log(airpass.df$passengers)~time + month, data = airpass.df )
anova(pass.fit1)
## Analysis of Variance Table
## 
## Response: log(airpass.df$passengers)
##            Df  Sum Sq Mean Sq  F value    Pr(>F)    
## time        1 25.1233 25.1233 7143.582 < 2.2e-16 ***
## month      11  2.2843  0.2077   59.047 < 2.2e-16 ***
## Residuals 131  0.4607  0.0035                       
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

p-value可以看成指示性的,因为不满足独立性假设。

  • 看残差图
plot( residuals(pass.fit1)~fitted(pass.fit1))
lines(lowess(fitted(pass.fit1), residuals(pass.fit1)))

二次方的?

  • 重新拟合
pass.fit2 = lm( log(airpass.df$passengers)~time + month + I(time^2), data = airpass.df )
anova(pass.fit2)
## Analysis of Variance Table
## 
## Response: log(airpass.df$passengers)
##            Df  Sum Sq Mean Sq   F value    Pr(>F)    
## time        1 25.1233 25.1233 10813.900 < 2.2e-16 ***
## month      11  2.2843  0.2077    89.386 < 2.2e-16 ***
## I(time^2)   1  0.1587  0.1587    68.307 1.409e-13 ***
## Residuals 130  0.3020  0.0023                        
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
plot(residuals(pass.fit2)~fitted(pass.fit2))
lines(lowess(fitted(pass.fit2), residuals(pass.fit2)))

plot(residuals(pass.fit2), type = "l") #link the points lower case L
lines(lowess(fitted(pass.fit2), residuals(pass.fit2)))

* 这回看起来好多了,相对恒定而且没有趋势了 * 但残差对时间的图表现出非独立性 4. 残差是否有自相关性?

acf(residuals(pass.fit2))

相关性和自相关性

线性相关系数\(\rho\),通常被当做相关性,是衡量两个随机变量\(X 和 Y 之间的线性关系,其中 0<\rho<1\)

  • \(如果 \rho = 1 \,则X和Y是完全正相关\)
  • \(如果 \rho = -1 则X和Y是完全负相关\)
  • \(如果 \rho = 0 \,则X和Y是完全不相关\)

  • 我们通过样本相关系数\(r\)来估算相关系数\(\rho\)
  • 在简单线性回归中样本的相关系数: \[r = \rm sign(\beta_1) \sqrt{R^2}\]

  • 自相关性是变量跟其自身的相关性,但不完全正确,变量Y跟其自身的相关性为1,因为他跟他自己是完全正相关的
  • 这里的自相关是延迟h的自相关,是观测值Y和那些h个时间单位以后的Y的值的相关性, T 为时间序列的长度 \[\rm Cor[Y_t, Y_{t+h}]\]
  • 例如:h=1时,我们指的相关性是存在于时间t和时间t+1的观测值的相关性
  • 这里的相关性是指 \(Y_t 相对于 Y_{t+h}\)的 其中: t=1,…, T-2.

ACF 图

  • ACF图是自变量h=0,1,…,T-2,的一个函数图,通常没到T-2就被截断了,因为数据太少了
  • 垂线的高度代表我们感兴趣的变量(通常为残差)在时间t和t+h的值的相关性
  • 第一条垂直线的高度恒为1,因为\(\rm Cor[Y_t, Y_{t}] = 1\)
  • 第二条垂直线的高度是\(\rm Cor[Y_t, Y_{t+1}]\), lag-1 自相关叫1阶自相关
  • 第三条垂直线的高度是\(\rm Cor[Y_t, Y_{t+2}]\), lag-2自相关又称为二阶自相关,以此类推

  • 蓝色的水平虚线是我们认为自相关性I显著的位置。大致来讲,如果垂直线超过虚线,则\(H_0: \rho_h = 0\) 的 P-value将会小于0.05。(由于多重比较,即使没有相关性,有一两个自相关值稍微超过虚线也是正常的

回到前面如何识别自相关性的第四步

  1. 残差是否有自相关性?
acf(residuals(pass.fit2))

* 这是明显的自相关模式,有很多垂直线超过了虚线。我们遇到了自相关问题。 * 要解决这个问题,我们明确建立一个一阶自相关模型,通过引进一个延迟响应变量作为解释变量到模型中来解决这个问题。 * 即:\(Y_1\)用来解释\(Y_2\)\(Y_2\)用来解释\(Y_3 \cdots\), 以此类推 * 这意味着我们不能使用\(Y_1\)作为相应变量,因为没有\(Y_0\) log(passengers)[-1] 去掉第一项 * 这类模型称为延迟响应模型

pass.fit3 = lm( log(passengers)[-1] ~ time[-1] + month[-1] + I(time[-1]^2) + log(passengers)[-144], data = airpass.df)
  • 我们总要加上延迟相应变量作为模型的最后一项。
anova(pass.fit3)
## Analysis of Variance Table
## 
## Response: log(passengers)[-1]
##                        Df  Sum Sq Mean Sq  F value    Pr(>F)    
## time[-1]                1 24.4515 24.4515 19260.14 < 2.2e-16 ***
## month[-1]              11  2.2733  0.2067   162.79 < 2.2e-16 ***
## I(time[-1]^2)           1  0.1617  0.1617   127.35 < 2.2e-16 ***
## log(passengers)[-144]   1  0.1362  0.1362   107.25 < 2.2e-16 ***
## Residuals             128  0.1625  0.0013                       
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
  • P-Value很小,有力证明延迟相应项应该包含在模型中(假设我们的模型是正确的)
acf(residuals(pass.fit3))

  • 至少有少量的最高点超过了虚线,这可能是因为多重比较导致的
  • 留意一阶自相关建模是如何消除其他显著的自相关的
    • 因为如果A和B紧密相关,而B和C紧密相关,则A和C也可能是紧密相关的。从而消除一阶自相关的时候可能同时消除了更高阶的关系。
    • 有两个延迟(12 和 16)也超过了虚线。这个acf图有效的做了21个假设测试,它给我们图像化展示了前21个自相关性是否显著(在5%水平上)。我们不应该感到奇怪去掉自相关性后还有如果1到两个是显著的。
  • 正式的通用延迟响应模型如下: \[Y_t = \beta_0 + \beta_1 t + \beta_2 X_{t1} + \cdots + \beta_{p+1} X_{tp} + \rho Y_{t-1}, \quad t = 2, \cdots , T\]
    • 通过包含一个延迟的响应变量到解释变量中,明确的调整一阶自相关。
    • 延迟相应变量的系数\(\rho\)可以当作残差自相关的天然估值。
    • 这个模型可以用来做预测但不适合作超过往前1步的预测。
  • 如果没有修正正自相关就拟合一个线性回归模型,则我们对残差标准差的估值将低于其实际值。
  • 结果,将会高估\(R^2\),还有回归系数的t-test以及anova F-test 也将变得无效。

为什么要关注自相关性

  • 信息的遗失:自相关的结果是阻碍我们了解总体中的个体差异的真实大小
  • 举例:

总结:

先说说GLM

  • 对于广义线性模型,这里讨论的主要就是当响应变量Y不是连续正态分布的时候。
  • 除了Poisson count data 是 log(E[Y]) 以外,其他无一例外都是log(odds) 或者称为 logit(p):
  • glm 建模的区别:
    • count data: family = poisson
    • binary data: family = binomial
    • count of success: cbind(success, fail), family = binomial
    • proportion: family = binomial, weight = n(trials)
    • 相同点都要看residual deviance vs degree of freedom
  • 如何通过把它们变换形式后拟合线性模型,求出它们的系数,之后通过逆变换回归分析它们自身,
  • 即:通过系数我们了解到当解释变量发生+1、×0.1、×2、×10(numerical explanatory),或分别处在 level 1, level 2 的时候,相应变量将发生什么相应的变化。
  • 由于Y进行了变换,有时候他们变换后的形式本身也具有统计意义,比如:
  1. count response:log-linear model
    • log(E[Y])
    • \(fac_1 * x2\) 直接解释吧,只是体现了glm 相对lm的优势,考虑了方差随均值的增加而增加,而且不用去掉本来有意义的outlier。
    • \(fac_1 * fac_2\) exponentiate coef 可以得到\(\hat{OR}\) 能解释的东西就多了,把2-by-2表格颠来倒去看就行。如果没有关联系,还能分列,分行独立进行对比解释。
    • binary response或者说binomial data或者说logistic regression,它的相应变量变为了logit(p)
    • 即 log(odds),所以我们只要exponentiate 那些系数就可以知道odds的变化,这样就行了,我们并不需要了解它的binary response的变化。
    • predict 了以后记得要逆变换才能得到p概率
    • \(\frac{exp(logit(p))}{1+exp(logit(p))} = p\)
  2. count of successes
    • 为什么成功数数据符合二项分布?因为他们是给定次伯努利实验的总和。
    • 它的fitting 比较特殊 有个cbind(success, fail)~x1*x2 — 在这种例子中,我们少了个体信息。还要假设给定行列,那么成功概率是确定的。每一次实验都是独立的。
    • 败血症例子中,因为可以死亡的人数是有上限的,所以即使响应变量是个count,我们也不能用泊松分布
    • 通过变换
  3. response is a proportion
    • S 型曲线在统计学语言里称为逻辑曲线,采用的方法称为逻辑回归
    • \(\beta_0 影响曲线的垂直高度,\beta_1 的量级(绝对值) 影响曲线的陡峭程度,越大越陡,\beta_1 < 0 的时候曲线是反S形。\)
    • log(odds)
    • 书本上predictCount()出来的可以直接解释为成功概率,不需要再exp了
  • 这里有几个有趣的地方:关联列表就是一个two-way anova
  • 它们的区别全部在响应变量的形式上:
  • 对于响应变量是:
    1. 各种组合下响应变量的值:比如来不来上课及段考过没过的学生的期考成绩
    2. 各种组合出现的次数:count to contingency table
    3. 各种组合下的成功与否:binary response,这个组合下成功为1不成功为0
    4. 各种组合成功的次数:count of successes 每种组合下成功的次数

考点

好吧,重头戏登场了:

MA & Excu Summary procedure

  • MA
    1. what & why a model was fitted & refitted
    2. all assumption was satisfied
    3. the final model is:
  • Excutive Summary
    1. problem of interest
    2. p-value
    3. estimation & prediction
    4. total variance explained & good or not for prediction

1. Null model was fit

2. simple linear model

  • MA
    1. 相应变量和解释变量都是数值quantitative variables. so a simple linear model was fitted

3. quadratic model

  • MA
    1. scatter plot suggested curvature
    2. residual plot(simple-linear.fit) showes curvature
  • Excu Summary
    1. For a one mark increase in assignment score, the increase in the expected exam score was greater as the the assignment score increased, e.g., there is little difference in expected exam scores for those getting 7 or 8 in assignment, much bigger for those getting 17 or 18.

4. Two-sample t-test for mean

  • MA:
    1. linear model was fitted
    2. with equal variance assumed in bot group)

5. Multitplicative model:

  • Multiplicative model 的特点:
    • 从散点图plot/trendscatter 看:Y值有指数级下降的趋势,伴随着离散度也下降
    • 从residual plot看,
      1. a log-linear model was fitted
      2. 前面几章的特点:求响应变量随解释变量变化的关系。
        • 注意:截距在这里没有意义,无须解释,因为截距已经不代表均值,而是当\(x_2, x_3 \cdots , x_n\)=0 的时候的值,而这些值在数据中通常都会不会为0.
        • 比如段考成绩为0的学生的期末考试成绩,是没有意义的。
      3. prediction interval 和 confidence interval的区别:predict的是预测样本中个体的值区间,confidence的是预测样本均值的区间,注意log-linear model只能预测median 而不是mean
      4. 对于相互影响的变量,我们通常需要知道的是它的调节作用,它调整了另一个变量的影响:It adjusts the effects of the test marks on exam marks.
    1. ANCOVA without interaction
    2. ANCOVA with interaction
    3. several explanatory variables.
      1. 通常将哪个变量被logged,哪个变量被去掉即可。
    4. Power law model. log-log model
    5. one-way ANOVA
      1. 一个分类变量,里面分类超过2中
    6. two-way ANOVA
      • MA
        • 是否parallel
        • 同A比B,同B比A 比大小趋势
        • 是否有 unusual high low
    7. poisson model
    8. poisson model
    9. logistic regression model
    10. time series
  • 残差的同分布,指的是残差图evocheck(fit)/plot(fit, which=1)
    • 图中对于每个相同x值上的多个y值(残差)在垂直方向上的分布是正态分布的,而不同x值对应的这个分布应该相同。所以应该都是分布在0周围,而且不应该有趋势:residuals (do not) appear to be randomly scattered around 0.

Exec Summary:

关于interpretation 的技巧:

  • 研究对象:
    • We want to quantify the relationship between … and …
    • We wish to investigate/check/quantify
    • We are interested in using…/checking whether…
  • 时态:
    1. 我们关注:we are interested …
    2. 图标显示:the plots showes …
    3. 我们估计(当时发生的情况): 用过去式
    4. 我们估计如果…会发生…: 可以用过去式也可以用一般现在时,will would
    5. between … and/to … 单位和百分比,可以写在后面一个数字后。
    6. 单位的单复数:表示量的多少需要复数,如5 grams, 1.8 meters. 如果是缩写则表示单位,不用复数:34.1 cm
  • 分段解释:
    1. A child has a higher epxected birth-weight the longer its gestiataion time - up to a 42 weeks - then it starts decreseing in size the longer it stays unborn.
    2. 这里的分段位置在42 weeks,之前正相关,之后负相关。
    3. 通过增加一个dummy variable,我们得到两个又相互影响的解释变量。
    4. dummy 变量影响另一个解释变量的效果。e.g. for non-regular attenders each additional mark of … the … increase somewhere between … and …, for students who attended regularly, this effect on the exam mark increased ()
  • \(R^2\) 解释的是数据的或或者是响应变量的variability 变异性
  • 快速反应:log(Y) 及 log-log的区别
    1. 对于:log-linear model
      • \(log(y_i) = \beta_0 + \beta_1 x_i\)
      • \(y = exp^{\beta_0} \times exp^{\beta_1 \cdot x}\)
      • \(for \, x_i^{'} = x_i+a\)
      • \(y_i^{'} = y_i \times exp^{\beta_1 \times a}\)
      • 所以x增加a,y就乘以 \(e^{\beta_1 \cdot a}\)
    2. 对于: log-log model
      • \(log(y_i) = \beta_0 + \beta_1 log(x_i)\)
      • \(y_i = exp^{\beta_0} \times exp^{\beta_1 \cdot log(x_i)}\)
      • \(y_i = exp^{\beta_0} \times x^{\beta_1}\)
      • \(for \, x_i^{'} = x_i \times a\)
      • \(y_i^{'} = y_i \times a^{\beta_1}\)
      • x乘以a,y就乘以 \(a^{\beta_1}\)
    3. 区别:
      • 注意在log-linear model中对比的x的增加量 为了让影响更加显著 \(\Delta_x= +1, +10, +100 \Rightarrow \Delta_y = \times e^{\beta_1 / 10\beta_1/100\beta_1}\)
      • 在log-log model中对比的是x变化的倍数\(\Delta_x= \times2, \times4, \times6 \Rightarrow \Delta_y = \times {(2/4/6)}^{\beta_1}\)
  • 常用句型:
    • seems to be/appears to be/does not appear to be
    • slight trend/variability increases with increaseing fitted values(y)
    • there is also …
    • roughly constant
    • more-or-less averaging around zero
    • come to a erroneous conclusion base on wishful thinking
  • 统计显著:
    • the P-value associated with the Attend variable is highly statistically significant.
  • Variation -差异 正常预期结果与观测结果之间的差额总量即为差异,variation表示变化,variance表示差异
  • 标准差
    • 简单来说,标准差是一组数据平均值分散程度的一种度量。一个较大的标准差,代表大部分数值和其平均值之间差异较大;一个较小的标准差,代表这些数值较接近平均值。

    \[\sigma = \sqrt{ \frac{1}{N} \sum_{i=0}^N (x_i-\mu)^2}\]

  • 方差
  • 标准误
    • 标准误差针对样本统计量而言,是某个样本统计量的标准差。当谈及标准误差时,一般须指明对应的样本统计量才有意义。以下以样本均值(样本均值是一种样本统计量)作为例子:
    • 例如, 样本均值是总体均值的无偏估计。但是,来自同一总量的不同样本可能有不同的均值。于是,假设可以从总体中随机选取无限的大小相同的样本,那每个样本都可以有一个样本均值。依此法可以到一个由无限多样本均值组成的总体,该总体的标准差即为标准误差。
    • \(SE_{\overline x} = \frac{s}{\sqrt n}\)
  • t值

  • iid:
    1. independence: investigate how we obtained the data.
    2. identically distributed: result in the variance of the residual being roughly constant(regardless of the fitted values) and the residuals more-or-less averaging around zero. (see with evocheck)
    3. Normality (see using normcheck)
  • 当要你comment一个plot的时候通常要你comment什么?
    1. yx1x2x3关系,甚至xx关系(相关性),
    2. 变换,曲线,转折点
    3. 异常点,极值点